* Figure: INTERVENTION STOVE PURCHASE RATES

use "${output}panel_r0_r1_r2.dta", clear // Load merged analysis panel
set seed 293847

********************************************************************************

// Keep households surveyed in each round
keep if sample_household == 1

// Reposition baseline value of household controls
foreach var of varlist household_size children_under_five {
	gen b_`var' = `var' if surveyround == 0
	bys hh_id (surveyround) : replace b_`var' = b_`var'[1] if mi(b_`var')
	}

// Merge with NGO activity raw data
merge m:1 gpname using "${data}chirag_activity_list.dta", ///
	keepusing(chirag_village health education forestry fodder spring youth_club agricultural herbs shg cooperative pirna_ws year_start)

// Sort uniquely for simulations
sort hh_id surveyround

// Keep relevant NGO strata variable
tab chirag_strata chirag_village, miss
drop chirag_village

// Number of simulations
local reps = 10000 // number of bootstrap repetitions

********************************************************************************

* COUNT OF NUMBER OF ACTIVITIES CONDUCTED BY THE NGO

* Generate count NGO-activity intensity count
egen ngo_activity_count = rowtotal(health education forestry fodder spring youth_club agricultural herbs shg cooperative pirna_ws), missing
lab var ngo_activity_count "Count of projects conducted by NGO (maximum possible = 11)"
replace ngo_activity_count = 0 if chirag_strata == 0
gen missing_ngo_activity_count = 0
replace missing_ngo_activity_count = 1 if mi(ngo_activity_count)
lab var missing_ngo_activity_count "=1 if no data available on range of activities"

* Bootstrap with randomly generated NGO activity counts for missing GPs
capture program drop ngo_count_boot_rep
program define ngo_count_boot_rep, rclass
	preserve
		gen ngo_activity_count_replace = runiformint(0, 7)
		replace ngo_activity_count = ngo_activity_count_replace if missing_ngo_activity_count == 1
		bsample, cluster(uniquegrp)
		reghdfe purchased_intervention_stove i.treatment##c.ngo_activity_count##c.ngo_activity_count b_* ///
			if surveyround == 1, absorb(districtcode) vce(cluster uniquegrp)
		return scalar beta_count_interaction = _b[1.treatment#c.ngo_activity_count]
		return scalar beta_count_interaction_sq = _b[1.treatment#c.ngo_activity_count#c.ngo_activity_count]
	restore
end

* Do the simulation and generate density plot for full sample
preserve

	qui simulate beta_ngo_count_boot_rep = r(beta_count_interaction) ///
		beta_ngo_count_sq_boot_rep = r(beta_count_interaction_sq), ///
		reps(`reps') nolegend seed(213344) : ngo_count_boot_rep

	lab var beta_ngo_count_boot_rep "{&beta}(TREATMENT{sub:j} × NGO{sub:j})"
	lab var beta_ngo_count_sq_boot_rep "{&beta}(TREATMENT{sub:j} × NGO{sup:2}{sub:j})"

	gen n = _n
	tempfile count_boot_rep
	save `count_boot_rep'

	foreach var of varlist beta_ngo_count_boot_rep beta_ngo_count_sq_boot_rep {

		_pctile `var', percentiles(2.5 5 95 97.5)
		
		graph twoway (kdensity `var', ///
			xtitle("`: variable label `var''") ytitle("") ///
			title("")) ///
			|| (kdensity `var', range(`r(r1)' `r(r4)') color(gs14%50) recast(area)) ///
			|| (kdensity `var', range(`r(r2)' `r(r3)') color(gs11%50) recast(area) ///
			legend(order(3 2) ///
				label(2 "95% C.I.") label(3 "90% C.I.") ///
				cols(2)) ///
			name(`var', replace) nodraw ///
			xline(0, lpattern(dash)))  

	}

	grc1leg beta_ngo_count_boot_rep ///
		beta_ngo_count_sq_boot_rep, ///
			leg(beta_ngo_count_boot_rep) ///
			title("(a) Count of NGO projects", size(small)) ///
			pos(3) cols(3) name(activities, replace)
	gr close

restore

********************************************************************************

* COUNT OF NUMBER OF YEARS OF NGO ACTIVITY

gen years_active = 2013 - year_start
replace years_active = 0 if chirag_strata == 0 
lab var years_active "Number of years NGO has been active in this GP (as of 2013)"
gen missing_years_active = 0
replace missing_years_active = 1 if mi(years_active)
lab var missing_years_active "=1 if no data available on starting year of NGO activities"

* Bootstrap with randomly generated NGO active years for missing GPs
capture program drop ngo_year_boot_rep
program define ngo_year_boot_rep, rclass
	preserve
		gen years_active_replace = runiformint(0, 25)
		replace years_active = years_active_replace if missing_years_active == 1
		bsample, cluster(uniquegrp)
		reghdfe purchased_intervention_stove i.treatment##c.years_active##c.years_active b_* ///
			if surveyround == 1, absorb(districtcode) vce(cluster uniquegrp)
		return scalar beta_year_interaction = _b[1.treatment#c.years_active]
		return scalar beta_year_interaction_sq = _b[1.treatment#c.years_active#c.years_active]
	restore
end

* Do the simulation and generate density plot for full sample
preserve	

	qui simulate beta_ngo_year_boot_rep = r(beta_year_interaction) ///
		beta_ngo_year_sq_boot_rep = r(beta_year_interaction_sq), ///
			reps(`reps') nolegend seed(348723) : ngo_year_boot_rep

	lab var beta_ngo_year_boot_rep "{&beta}(TREATMENT{sub:j} × NGO{sub:j})"
	lab var beta_ngo_year_sq_boot_rep "{&beta}(TREATMENT{sub:j} × NGO{sup:2}{sub:j})"

	gen n = _n
	tempfile year_boot_rep
	save `year_boot_rep'

	foreach var of varlist beta_ngo_year_boot_rep beta_ngo_year_sq_boot_rep  {

		_pctile `var', percentiles(2.5 5 95 97.5)

		graph twoway (kdensity `var', ///
			xtitle("`: variable label `var''") ytitle("") ///
			title("")) ///
			|| (kdensity `var', range(`r(r1)' `r(r4)') color(gs14%50) recast(area)) ///
			|| (kdensity `var', range(`r(r2)' `r(r3)') color(gs11%50) recast(area) ///
			legend(order(3 2) ///
				label(2 "95% CI") ///
				label(3 "90% CI") ///
				cols(2)) ///
			name(`var', replace) nodraw ///
			xline(0, lpattern(dash)))		
	}

	grc1leg beta_ngo_year_boot_rep ///
		beta_ngo_year_sq_boot_rep, ///
			leg(beta_ngo_year_boot_rep) ///
			pos(3) cols(2) name(years, replace) ///
			title("(b) Years of NGO activity", size(small))
	gr close

restore

********************************************************************************

grc1leg activities years, ///
	legendfrom(activities) cols(1) pos(6) iscale(*1.2) ///
	l1("Density", size(small)) name(allvillages_rep_quad, replace)

gr disp allvillages_rep_quad, ysize(12) xsize(10)

graph export "${results}figure_ngo_redefinition_quadratic.png", ///
	width(5000) replace

gr close

********************************************************************************

use `count_boot_rep', clear
merge 1:1 n using `year_boot_rep'

* END
